Embedding methods for large-scale surface calculations 
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One of the goals in the development of large scale electronic structure methods is to perform 
calculations explicitly for a localised region of a system, while still taking into account the rest 
of the system outside of this region. An example of this in surface physics would be to embed 
an adsorbate and a few surface atoms into an extended substrate, hence considerably reducing 
computational costs. Here we apply the constrained electron density method of embedding a Kohn- 
Sham system in a substrate system (first described by P. Cortona[l|] and T.A. WesolowskiQ]), within 
a plane-wave basis and pseudopotential framework. This approach divides the charge density of the 
system into substrate and embedded charge densities, the sum of which is the charge density of the 
actual system of interest. Two test cases are considered. First we construct fee bulk aluminium by 
embedding one cubic lattice of atoms within another. Second, we examine a model surface/adsorbate 
system of aluminium on aluminium and compare with full Kohn-Sham results. 

PACS numbers: 71.15.Mb,71.15.Ap,71.15.-m 
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I. DFT EMBEDDING 

Density Functional Theory (DFT) is one of the most 
powerful tools for the ab initio calculation of the physi- 
cal and chemical properties of materials, being both effi- 
cient and accurate. Many implementations of DFT exist, 
which differ in the approach taken to approximating the 
unknown density functional that describes the contribu- 
tion to the energy of the electrons that is not due to 
the external potential. The basic problem of DFT is to 
minimise the functional 



E[p] = T s \p] + J[p] + E xc [p\ + J V ext {r)p{r)d 3 



(1) 



where T s [p] is the non-interacting kinetic energy func- 
tional, J[p] is the Hartree energy, E xc [p] is the 'exchange- 
correlation' energy which takes into all non-classical 
electron-electron interactions, and V ex t(r) is the exter- 
nal potential of the system of interest. To obtain the 
ground state energy this must be minimised with respect 
to variations in p(r), subject a constrained number of 
electrons. The Hartree and external potential energies 
can be simply evaluated for a trial charge density, and 
accurate approximations are available for E xc [p], but no 
explicit form is known for the non-interacting kinetic en- 
ergy. The charge density for the interacting electrons 
can be identified with the charge density of a 'reference' 
non-interacting electron gas, and this reference system 
solved self-consistently to yield the energy and charge 
density for the interacting system - this is the Kohn-Sham 
method. Due to orthogonalisation requirements, solving 
for this reference system scales as 0(N 3 ) where N is the 
number of electrons which limits the size of system that 
can be considered Q. 



There are O(N) methods available that take advan- 
tage of the 'nearsightedness' of the density matrix, but 
at present these are only applicable where the number of 
basis functions per atom is small, such as tight-binding 
or atomic-orbital calculations. Another approach is to 
employ approximate kinetic energy functionals, and min- 
imise the functional directly. Calculations of this sort are 
cheap and quick, but the approximations available in the 
literature are generally not sufficient for structural opti- 
misation, let alone chemical accuracy. 

Another approach is embedding. In many cases (eg 
an adsorbate on a surface) we can divide the system into 
two regions of space, region / (eg the adsorbate and a few 
surface atoms) and region 77 (eg the rest of the surface). 
Region II is largely the same as that for a more simple 
system that may easily be solved for, whereas region / 
is where all of the interesting physics occurs. It would 
obviously be computationally advantageous to solve for 
region II first, and then solve for region / taking into 
account the influence of region II in some way. This 
'embedding' approach has received a great deal of atten- 
tion, and a large number of methods have been presented 
in the literature Many methods approach this as 

a boundary value problem, at the wavefunction level, but 
we examine a different approach which starts at the more 
flexible DFT level, and was first presented by Cortona[lj 
and WesolowskiQ- 

We start with Eq. (1) written as 

E[p] = T s [ Pl }+Tr dd [Pi,P2}+T s [p 2 ] 

+ J[p] + E xc [p] + [ V ext {r)p(r)d 3 r (2) 



where p(r) — pi(r) +pa(r), Pi(r) is the charge density of 
the embedded system that we shall be varying, and p 2 ('r) 
is the substrate charge density that is kept constant. This 

, pz] , as 



defines the non- additive kinetic energy, T™ add [pi, 
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Tmadd r rt 
s [Pi: 



P%] = T s [pi + p 2 ] - T s [pi] - T s [p 2 



(3) 
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Minimising Eq. (2) with respect to variations in pi(r) 
leads to the Euler-Lagrange equation 



5T s [ Pl 
Spi 
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[pi, Pi] 



8 pi 



+ V K s[p;r]= l i (4) 



where p is the chemical potential, Vks [P'i r ] is the usual 
Kohn-Sham potential associated with density p(r), and a 
new 'embedding potential' term is present. In the same 
manner as for the Kohn-Sham case, this leads to the 
pi(r) being the solution of the 'Kohn-Sham' equations 
associated with Eq. (4) at self consistency, but with an 
effective potential given by 



V e "[p;r] = V KS [p\v] 
= V K s[p\ r] 
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(5) 



As applied previously, this method has employed basis 
function localised to their host atoms, which essentially 
constrains the charge density to a localised region and 
to a particular functional form. These previous appli- 
cations have been to weakly interacting and insulating 
systems, whereas here we explore the possibility of ex- 
amining strongly interacting metallic systems using an 
unbiased plane wave basis. 



II. APPROXIMATE KINETIC ENERGY 
FUNCTION ALS 

As given above we have only re-expressed the original 
problem in a slightly different form. In order for this ap- 
proach to be useful, it must be possible to express the 
non-additive kinetic-energy accurately. Of course we do 
not know an analytic form for this interaction energy, 
but it is expected to be small, and zero for no overlap 
between the embedded and substrate systems. This sug- 
gests the use of the approximate kinetic energy function- 
als available in the literature to construct T™ add and its 
functional derivative. 

Two forms of approximate kinetic energy functional 
are considered. First a semi-local form, incorporating an 
'enhancement factor' analogous to the Generalised Gradi- 
ent Approximation (GGA) of exchange-correlation func- 
tionals. This takes the form 



l [p\ = §(3^)1 



^F(t)d 3 



where 
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(6) 



(7) 



and F(t) is the enhancement factor. The Thomas-Fermi 
functional is a special case of this, where F(t) — 1, the 
I s * order gradient expansion corresponds to F(t) = 1+at 
with a an appropriate constant, and the von Weizacker 



approximation corresponds to F(t) — bt with b again an 
appropriate constant. 

The function F(t) is generally chosen to provide appro- 
priate limiting behaviour for the functional, and parame- 
ters are often chosen to fit data, theoretical or experimen- 
tal. The functional derivative of this can be expressed 
most concisely as 
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(8) 



which differs from the expression obtained by a direct ap- 
plication of the usual formulae Q. Severe aliasing prob- 
lems arise if the standard form is applied, as occurs for 
the exchange-correlation potential but this alternative 
analytic form greatly reduces these numerical difficulties. 

The above functional has many deficiencies, the pri- 
mary one being of only limited non-locality. Explicitly 
non-local functional have been investigated in the liter- 
ature, with the introduction of an analytic form that 
integrates over the contributions from the charge den- 
sity at each pair of points in real space @. This can be 
further generalised by adding a third order term that 
integrates the contributions from triplets of points in 
space, and higher order terms. Unfortunately these are 
generally extremely computationally expensive to evalu- 
ate. One exception to this is a form proposed by W ang 
and Teter[l(|, Perrotp4| and Smargiassi and Madden [l2f 
where the non-local term is expressed as a convolution 
integral which can be evaluated efficiently in reciprocal 
space. This approximate functional is given by 
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pa V 2 p?d 3 r 
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+ J p Q (r)w Q (r-r> Q (r')d 3 rdV (9) 

where the first term is the Thomas Fermi contribution, 
the second the von Weizacker contribution, and the final 
term a non-local contribution. The parameter a is arbi- 
trary, and w a (r — r') is chosen such that the functional 
has the correct linear response for a homogeneous non- 
interacting electron gas with the same average charge 
density as the charge density of interest. The functional 
derivative of this approximation is given by 
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(37T)3p3 - p- 



+2ap a ~ 1 {r) I p a (r')w a (r-r')d 3 r' 



(10) 



Our aim is to assess the usefulness of these approxima- 
tions for carrying out the embedding procedure described 
in the previous section. 
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Functional E[p]/eV AE/eV R/% 
Tpwsg -59.68 -1.35 7.2 

T nlcc _ 5g41 _ 00 g 12 

2 

Kohn-Sham -58.33 



TABLE I: Total energies per atom, and errors in energies and 
charge density. 



III. RESULTS 

To investigate this partially frozen density approach we 
examine the first test case, fee aluminium with a conven- 
tional 4 atom cubic unit cell. The 3 face-centred atoms 
are taken to be the substrate system, and this structure 
is solved for first to provide P2(r) and T s [p 2 ]- A stan- 
dard plane- wave method [l3| is used, with a lattice con- 
stant of do = 4.05A, a plane-wave cut-off of 200 eV, 35k 
points in the irreducible wedge, and the Goodwin-Needs- 
Heine local pseudopotential,14j. Exchange-correlation is 
described by the LDA. 

Once this substrate is constructed the embedded 
Kohn-Sham calculation is carried out as a standard 
plane-wave calculation, but with the trial potential given 
by Eq. (5) and the total energy given by Eq. (2). These 
additional terms describing the substrate, and the cal- 
culation only involves the three electrons introduced by 
the embedded atom at the corner of the unit cell; the 
9 substrate electrons are taken into account entirely by 
their charge density. Parameters of the calculation are 
chosen to be the same as for the substrate calculation. 
It should be made clear that for the substrate calcula- 
tion we are solving for the lattice of face centred atoms 
and their accompanying electrons, but for the embedded 
calculation we are solving for the entire fee system, but 
only the electrons associated with the embedded (corner) 
atom are provided with a Kohn-Sham representation. 

We discuss results for an enhancement factor approx- 
imation (Eq. (6)) with the Perdew and Wang '86 en- 
hancement factor [l5j (Tpy/m)- 

F(s) = (1 + 1.296s 2 + 14s 4 + 0.2s 6 )ts 



and the non-local linear-response corrected functional, 
Eq. (9) with a = \ (Tf oc ). Although calculations have 

been carried out for a large number of different enhance- 
ment factor functionals there was no significant differ- 
ence in the results and any small differences do not affect 
our conclusions. Table 1 shows the errors in the total 
energies, and the charge density expressed as the mean 
absolute deviation as a percentage of the mean density. 
It is immediately apparent that the non-local functional 
provides the superior approximation. 

From this we conclude that by applying the non-local 
functional an accurate total energy and charge density 



can be obtained by this DFT embedding procedure for 
an essentially metallic and strongly interacting system. 
Semi-local enhancement factor functional do not result 
in a useful accuracy. 

Next we move onto the second test model sur- 

face/adsorbate system, the type of system we hope to 
apply the method to in the future. This system also 
places far more demand on the method since the charge 
densities are considerably more inhomogeneous. We con- 
sider a \J 7 l x \[2 super-cell with one 'adsorbate' per cell, 
a 3-layer (100) slab and 5 equivalent vacuum layers. The 
adsorbate is centred on a four-fold hollow site. The sub- 
strate (p2(r)) is chosen to be the lower two layers, and 
embedding calculations were performed with three em- 
bedded atoms making up the upper surface layer and 
the adsorbate. All calculations were performed with 10k 
points in the irreducible wedge. We chose to examine the 
potential energy curves produced by varying the adsor- 
bate/surface, and compare these with the equivalent full 
Kohn-Sham results. We also compare the embedding re- 
sults with Kohn-Sham results for a 1-layer (100) slab to 
discover whether the embedding method can reproduce 
the interaction between the adsorbate + top layer and 
lower two layers. 
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FIG. 1: Adsorption energy as a function of surface/adsorbate 
distance vertically above a four-fold hollow. Dashed line is 3- 
layer + adsorbate Kohn-Sham result. Dotted line is 1-layer + 
adsorbate Kohn-Sham result. Solid line is results for embed- 
ding calculation. Zero energy is sum of surface and isolated 
adsorbate total energies. 

Figure 1 shows the potential energy curves for the em- 
bedded calculation with the non-local functional, for a 
full Kohn-Sham calculation and for a Kohn-Sham calcu- 
lation incorporating only one layer of the surface. Results 
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are promising, with the embedding results agreeing well 
with the full 3-layer Kohn-Sham calculation. 

Equilibrium results for the embedded system are a sur- 
face/adsorbate distance of 1.82 A and —2.07 eV, com- 
pared to Kohn-Sham results of 1.77 A and -2.04 eV - 
an error of 0.05 A and 0.03 eV respectively. This differs 
greatly from the 1 layer/ adsorbate system (1.48 A and 
—2.26 eV). No results are presented for semi-local en- 
hancement factor functionals, since these differ negligibly 
from the Kohn-Sham 1 layer/adsorbatc layer results, in- 
dicating that these functionals cannot take into account 
the coupling between the substrate and adsorbate + top 



layer for this system. 

We conclude that the DFT embedding approach, with 
a non-local functional, can provide accurate results for 
two systems with very different charge densities. Both 
of these systems are metallic and strongly interacting, 
so this approach may be useful for investigating a wide 
range of large systems. 

Future work will consist of the application of this 
method to further systems, both for further validation 
of the accuracy and the study of large adsorbate/surface 
systems. 
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